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Abstract. This paper presents a new approach for analysing structural properties of 
time series from complex systems. Starting from the concept of recurrences in phase 
space, the recurrence matrix of a time series is interpreted as the adjacency matrix of 
an associated complex network which links different points in time if the evolution of 
the considered states is very similar. A critical comparison of these recurrence networks 
with similar existing techniques is presented, revealing strong conceptual benefits of 
the new approach which can be considered as a unifying framework for transforming 
time series into complex networks that also includes other methods as special cases. 

It is demonstrated that there are fundamental relationships between the topological 
properties of recurrence networks and the statistical properties of the phase space 
density of the underlying dynamical system. Hence, the network description yields 
new quantitative characteristics of the dynamical complexity of a time series, which 
substantially complement existing measures of recurrence quantification analysis. 
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1. Introduction 

Since the early stages of quantitative nonlinear sciences, numerous conceptual 
approaches have been introduced for studying the characteristic features of dynamical 
systems based on observational time series [1-4]. Popular methods that are increasingly 
used in a variety of applications (see, for example, [5]) include (among others) Lyapunov 
exponents, fractal dimensions, symbolic discretisation, and measures of complexity such 
as entropies and quantities derived from them. All these techniques have in common 
that they quantify certain dynamically invariant phase space properties of the considered 
system based on temporally discretised realisations of individual trajectories. 

As a particular concept the basic ideas of which originated in the pioneering work 
of Poincare in the late 19th century [6], the quantification of recurrence properties in 
phase space has recently attracted considerable interest [7]. One particular reason for 
this is that these recurrences can be easily visualised (and subsequently quantified in a 
natural way) by means of so-called recurrence plots obtained from a single trajectory of 
the dynamical system under study [8,9]. When observing this trajectory as a scalar time 
series x{t) {t = 1, . . . , N), one may use a suitable m-dimensional time delay embedding 
of x{t) with delay r [11], x.^"^\t) = {x{t), x{t + r), . . . , x{t + (m — l)r)), for obtaining a 
recurrence plot as a graphical representation of the binary recurrence matrix 



where G(-) is the Heaviside function, || ■ || denotes a suitable norm in the considered 
phase space, and e is a threshold distance that should be reasonably smaller than the 
attractor diameter [9, 10]. To simplify our notation, we have used the abbreviation 
Xj = x^^\t = ti) (with ti being the point in time associated with the i-th observation 
recorded in the time serie^ wherever appropriate. 

Experimental time series often yield a recurrence plot displaying complex structures, 
in particular, with different properties of the non-interrupted diagonal and vertical 
structures ("lines"). A variety of statistical characteristics of the length distributions 
of these lines (such as maximum, mean, or Shannon entropy) can be used for 
defining additional quantitative measures that characterise different aspects of dynamic 
complexity of the studied system in more detail. This conceptual framework is known 
as recurrence quantification analysis (RQA) [12-14] and is nowadays frequently applied 
to a variety of real-world applications of time series analysis in various fields of 
research [15]. However, most of these RQA measures are sensitive to the choice of 
embedding parameters, which are found to sometimes induce spurious correlations in a 
recurrence plot [16]. 

Recent studies have revealed that the fundamental invariant properties of a 
dynamical system (i.e., its correlation dimension D2 and correlation entropy K2) are 
conserved in the recurrence matrix [17]. Furthermore, it is found that the estimation of 

I Note that unlike many other methods of time series analysis, the concept of recurrence plots does 
not require observations that are equally spaced in time. 




(1) 
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these invariants is independent of the particular embedding parameters. The recurrence 
plots preserve all the topologically relevant phase space information of the system, such 
that one can completely reconstruct a time series from its recurrence matrix (modulo 
some rescaling of its probability distribution function) [18, 19]. 

A further appealing paradigm for analysing structural features of complex systems 
is based on their representation as complex networks of passive or active (i.e., mutually 
interacting) subsystems. For this purpose, classical graph theory has been systematically 
extended by a large variety of different statistical descriptors of the topological features 
of such networks on local, intermediate, and global scales [20-22]. These measures 
have been successfully applied for studying real-world networks in various scientific 
disciplines, including the structural properties of infrastructures [23], biological [24], 
ecological [25], and climate networks [26], to give some prominent examples. The 
corresponding results have triggered substantial progress in our understanding of the 
interplay between structure and dynamics of such complex networks [27-29]. 

The great success of network theory in various fields of research has recently 
motivated first attempts to generalise this concept for a direct application to time 
series [30-39]. However, a substantial number of the recently suggested techniques 
have certain conceptual limitations, which make them suitable only for dealing with 
distinct types of problems. As an alternative that may provide a unifying conceptual 
and practical framework for nonlinear time series analysis using complex networks, we 
reconsider the concept of recurrences in phase space for defining complex network 
structures directly based on time series. For this purpose, it is straightforward to 
interpret the recurrence matrix R(e) as the adjacency matrix A(e) of an unweighted and 
undirected complex network, which we suggest to call the recurrence network associated 
with a given time series. To be more specific, the associated adjacency matrix is given 

by 

A,j(e)^R,,j(e)-5i,j, (2) 

where 5ij is the Kronecker delta introduced here in order to avoid artificial self-loops. 
A corresponding conceptual idea has recently been independently suggested by different 
authors [32,33,36,37,39], but not yet systematically studied. In this work, we however 
aim to give a rigorous derivation and detailed interpretation for a variety of quantitative 
characteristics of recurrence networks. It shall be noted that a generalisation to weighted 
networks (as partially studied in [30,31]) is straightforward if the recurrence matrix 
is replaced by the associated distance matrix between pairs of states. In any case, 
recurrence networks referring to the mutual phase space distances of observational points 
on a single trajectory are spatial networks, i.e. fully embedded into an m-dimensional 
space, which has important implications for their specific topological features. We will 
raise this point in detail within this paper. 

The consideration of recurrence plots as graphical representations of complex 
networks allows a reinterpretation of many network-theoretic measures in terms of 
characteristic phase space properties of a dynamical system. According to the ergodicity 
hypothesis, we suppose that one may gain full information about these properties by 
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either ensembles of trajectories, or sufficiently long observations of a single trajectory. 
Following this line of ideas, we may approximate the (usually unknown) invariant density 
p(x) (which is related with the associated invariant measure n by dn = p{x.)dx) of 
the studied system by some empirical estimate p^^\x) obtained from a time series, 
where e defines the level of coarse-graining of phase space involved in this procedure. 
Transforming the time series into a recurrence network then allows to quantitatively 
characterise the higher-order statistical properties of the invariant density p(x) by means 
of complementary methods, i.e., network-theoretic measures. 

According to the above argumentation, quantitative descriptors of the topological 
features of recurrence networks can be considered as novel measures within the 
framework of RQA. Our technique therefore exhibits additional deep insights into 
the phase space properties of dynamical systems directly related to their complex 
dynamics. Additionally, we will emphasise that we may also take seriously the duality 
of adjacency matrices of complex networks on the one hand, and recurrence matrices 
of dynamical systems on the other hand, which would also allow transferring concepts 
from dynamical systems theory (given that the corresponding recurrence plot based 
estimates are invariant under temporal reordering) to complex network theory. In this 
work, however, we will concentrate on a detailed discussion of how phase space properties 
can be further quantified in terms of network theory. 

The remainder of this paper is organised as follows: Section |2] presents a critical 
review of existing approaches for extracting complex networks from time series, 
including a comprehensive discussion of their potentials and potential problems (with 
a special emphasis on how to interpret the resulting networks' topological properties). 
The concept of recurrence networks as a natural alternative is further discussed in 
Section [31 In particular, it is demonstrated that many network-theoretic measures 
yield sophisticated quantitative characteristics corresponding to certain phase space 
properties of a dynamical system that have not yet been explicitly studied in terms 
of other dynamical invariants or measures of complexity based on RQA. In order to 
support our theoretical considerations. Section Hlprovides some examples of how different 
network-theoretic measures reveal certain phase space properties of various dynamical 
systems. Finally, we summarise our main results and outline some future directions of 
further methodological developments based on our proposed technique. 

2. Approaches for transforming time series into complex networks - A 
comparative review 

In this section, a review and classification of existing approaches for studying the 
properties of time series by means of complex networks approaches is presented (see 
Tab. [1]). In particular, strengths and possible limitations of the existing approaches will 
be briefiy discussed. 
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Method 



Vertex 



Edge 



Coarse-graining (2005) [44,49] 
Cycle networks (2006) [30] 
Correlation method (2008) [34] 
Visibility graph (2008) [35] 
Neighbourhood network (2008) [32,33] 
Recurrence network (2008) [36, 39] 



Discrete state s{t) 



Scalar state x{t) 



State x(t) 
State x(t) 



Cycle 
State x(t) 



Equality of states 
Correlation between cycles 
Correlation between state vectors 
Mutual visibility of states 
Recurrence of states (Mass) 
Recurrence of states (Volume) 



Table 1. Summary of the definitions of vertices and the criteria for the existence 
of edges in existing complex network approaches to time series analysis (given in 
chronological order). 



2.1. Coarse- graining of phase space 

The simplest possible method for transforming a time series into a complex network 
representation is coarse-graining its range into a suitable set of classes and considering 
the transition probabilities between these classes in terms of a weighted network. In 
general, the underlying concept of symbolic dynamics [40] allows characterising the 
properties of dynamical systems based on a partitioning of its phase space, yielding 
a transformation of every possible trajectory into an inifinite sequence of abstract 
symbols. Formally, the application of resulting methods of symbolic time series 
analysis (such as mutual information or entropic quantities) requires the existence of a 
generating partition which corresponds to a unique assignment of symbolic sequences 
(i.e., sequences of class identifiers) to every trajectory of the system. Note that this 
prerequisite is usually violated in real-world applications due to the presence of noise, 
however, even in the ideal noise-free case generating partitions do either not exist or 
can hardly be estimated (see [41] and references therein). Nevertheless, applications of 
symbolic time series analysis have recently attracted considerable interest in numerous 
applications [42,43]. 

Gao et al. [44-47] used a specific coarse-graining of the phase space for studying 
transitions between traffic states in different cellular automaton models. They however 
restricted their attention to the consideration of degree distributions. Recently, their 
idea was generalised to weighted networks by Zheng and Gao [48]. In a similar way, 
Li et al. [49, 50] used a corresponding approach for a coarse-grained analysis of stock 
exchange time series. In particular, their research mainly focussed on the identification 
of vertices with the highest relevance for information transfer, which have been quantified 
in terms of betweenness centrality and inverse participation ratio of the individual 
network vertices. 

The main disadvantage of the coarse-graining approach is that it may yield 
a significant loss of information on small amplitude variations. In particular, two 
observations with even very similar values may not be considered to belong to the 
same class if they arc just separated by a class boundary. This might infiuence the 
quantitative features of a corresponding network, as it is not exclusively determined by 
the widths of the individual classes, but also their specific definition. In this respect. 
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the recurrence network approach introduced in this work is more objective as it only 
depends on a single parameter e. Note, however, that coarse-graining might be a valid 
approach in case of noisy real-world time series, where extraction of dynamically relevant 
information hidden by the noise can be supported by grouping the data. 

2.2. Cycle networks 

In 2006, Zhang and Small [30-32] suggested to study the topological features of pseudo- 
periodic time series (representing, for example, the dynamics of chaotic oscillators like 
the Lorenz or Rossler systems) by means of complex networks. For this purpose, 
individual cycles (defined by minima or maxima of the studied time series) have been 
considered as vertices of a network, and the connectivity of pairs of vertices has been 
estabhshed by considering a generalisation of the correlation coefficient to cycles of 
possibly different length or, alternatively, their phase space distance. 

A potential point of criticism to this method is that the definition of a cycle is 
not necessarily straightforward in complex oscillatory systems. In [30-32], the authors 
mainly considered nonlinear oscillators in their phase- coherent regimes, however, it is not 
clear how a cycle could be defined for non-phase-coherent oscillations, for example, in the 
Funnel regime of the Rossler system. The same problem arises for systems with multiple 
time scales which are hence hard to treat this way. Furthermore, it is not intuitive how 
one can interpret correlations of cycles, since the values of the corresponding measures 
are not exclusively determined by the proximity of the corresponding parts of the 
trajectory in phase space, but depend also on the "lengths" of the cycles (in terms 
of the number of states) that may vary due to the discrete sampling. Hence, one may 
find rather low correlations between two cycles although the two parts of the trajectory 
are actually close to each other. 

2. 3. Correlation networks of embedded state vectors 

A generalisation of the technique used by Zhang and Small that can also be applied 
to time series without obvious oscillatory components has been suggested by Yang and 
Yang [34] using a simple embedding of an arbitrary time series. In their formalism, 
individual state vectors x^"^-* in the m-dimensional phase space of the embedded 
variables are considered as vertices, from which a Pearson correlation coefficient Vij = 
r(Xj^"*\ x^-"*^) can be easily computed. If rij is larger than a given threshold, vertices 
i and j are considered to be connected. The approach of Yang and Yang has recently 
been reinvented by Gao and Jin [39, 51] in terms of so-called fluid-dynamic complex 
networks (FDCN) that have been successfully applied for characterising the nonlinear 
dynamics of conductance fiuctuating signals in a gas-liquid two-phase fiow. 

One potential conceptual problem of this particular technique is that the 
consideration of correlation coefficients between two phase space vectors usually requires 
a sufficiently large embedding dimension m for a proper estimation with low uncertainty 
(more specifically, the standard error of the correlation coefficient is approximately 



Recurrence networks - A novel paradigm for nonlinear time series analysis 



7 



proportional to 1/y/m — 1). Hence, local information about the short-term dynamics 
captured in a time series might get lost when following this approach. Even more, since 
embedding is known to induce spurious correlations to a system under study, the results 
of the correlation method of network construction may suffer from related effects. 

With respect to the interpretation of the resulting network patterns, one has to 
note that for vertices corresponding to mutually overlapping time series segments, 
the consideration of correlation coefficients, as applied in both papers cited above, 
corresponds to studying the local auto-correlation function of the signal. Hence, the 
presence of edges between these vertices is exclusively determined by linear correlations 
within the signals. In principle, we might think of replacing the correlation coefficient 
by other measures of interrelationships such as the mutual information, that are 
also sensitive to general statistical dependences [26,52,53]; however, the appropriate 
estimation of such nonlinear quantities would require an even considerably larger amount 
of data, i.e., a very large embedding dimension m. 

Finally, when studying time series with pronounced cycles (like trajectories of the 
Rossler or Lorenz systems) , there may then be different cycles included in one embedding 
vector (depending on the sampling rate), which casts additional doubts with respect to 
the direct interpretability of the resulting network properties. 

2.4- Visibility graphs 

An alternative to the latter two threshold-based concepts has been suggested by Lacasa 
et al. [35] in terms of the so-called visibility graph. In this formalism, individual 
observations are considered as vertices, and edges are introduced whenever a partial 
convexity constraint is fulfilled, i.e. x{ta) and x{tb) are connected if for all states x{tc) 

with ta<tc<tb, 



holds. Visibility graphs have been used to study the behaviour of certain fractal as 
well as multifractal stochastic processes [54,55], energy dissipation in three-dimensional 
turbulence [56], and the nonlinear properties of exchange rate time series [57]. 

The network corresponding to a visibility graph is easily established and allows to 
distinguish between different types of systems. However, there is no straightforward 
interpretation of the convexity constraint in terms of phase space properties of the 
considered system. Moreover, the application of this approach is restricted to univariate 
time series, while at least the approach by Yang and Yang could in principle be easily 
generalised to multivariate time series. 

2. 5. Complex networks based on neighbourhood relations in phase space 




(3) 



As it has already been mentioned, the transformation of time series in terms of 
neighbourhood relationships has already been discussed by different authors. In 
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Recurrence network 



Phase space 



Vertex 

Edge 

Path 



State x(t) 
Recurrence of states 
Overlapping sequence of e-balls 



Table 2. Relationships between recurrence network entities and corresponding 
geometrical objects and their properties in phase space. 



particular, there are two possible approaches, that can directly be related to slightly 
different definitions of recurrence plots [9]: 

On the one hand, a neighbourhood can be defined by a fixed number of nearest 
neighbours of a single observation, i.e., a constant "mass" of the considered environments 
[32,33]. We refer to this method as a neighbourhood network in phase space. This 
setting implies that the degrees of all vertices in the network are kept fixed at the 
same value. Hence, information about the local geometry of the phase space, which is 
mainly determined by the invariant density p(x), cannot be directly obtained by most 
traditional complex network measures (see Sec. aS, Note that the adjacency matnx of 
a neighbourhood network is in general not symmetric, i.e., the fact that a vertex j is 
among the k nearest neighbours of a vertex i does not imply that i is also among the k 
nearest neighbours of j. Hence, neighbourhood networks can be formally considered as 
(partially) directed networks. 

On the other hand, one may define the neighbourhood of a single point in phase 
space (represented by a certain observation) by a fixed phase space distance, i.e., 
considering a constant "volume" [36,39]. This approach has the advantage that the 
degree centrality k^ gives direct information about the local phase space density (see 
Section [3]). Gao and Jin [39] termed a corresponding approach as fluid-structure complex 
networks (FSCN) and used it for analysing gas-liquid two-phase flow and the Lorenz 
system as a toy model in terms of link density. In addition, they related their results to 
the presence of unstable periodic orbits in a dynamical system. We will come back to 
this point in Section 14.31 

The consideration of neighbourhood relationships within a fixed phase space 
volume corresponds to the standard definition of a recurrence plot as mentioned in 
the introductory section. Hence, the resulting networks will be referred to as recurrence 
networks in the following. In particular, all arguments provided in the remainder of this 
paper for recurrence networks are based on the idea of a fixed volume of the considered 
neighbourhoods rather than a fixed mass and may not be directly generalised to the 
other case. 
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3. Quantitative assessment of recurrence networks 

Many of the already existing methods for transforming time series into complex networks 
that have been discussed in Section [2] suffer (among other problems) from the fact that 
there is no direct link between the local properties of the considered time series and the 
topology of the resulting complex networks. In particular, the concepts used for defining 
both vertices and edges of the networks, which differ across the various techniques, are 
in some cases rather artificial from a dynamical systems point of view (Table [I]). 

In contrast to the other recently suggested approaches, the identification of a 
recurrence matrix with the adjacency matrix of a complex network is a straightforward 
and natural idea that conserves many local properties of the considered time series. In 
particular, individual values of the respective observable can be directly considered as 
vertices of the recurrence network (similar to the visibility graph concept), while the 
existence of an edge serves as an indicator of a recurrence, i.e., pairs of states whose 
values do not differ by more than a small value e in terms of a suitable norm in phase 
space. 

It should be noted that the recurrence networks approach followed in this work 
is not the only concept that combines basic ideas of recurrence plots and complex 
networks. Apart from the neighbourhood networks originated in the idea of a fixed 
recurrence rate (i.e., a fixed mass of the considered neighbourhoods), the idea of 
considering a threshold value to the proximity of two vertices can also be found in other 
previously suggested methods. In particular, the coarse-graining approach is equivalent 
to considering recurrence plots of discrete-valued observables with a threshold of e = 0. 
Moreover, the correlation method of Yang and Yang [34] (see Sec. 12. 3p can also be 
considered as being based on a recurrence plot where the usual metric distance has been 
replaced by the correlation distance [58] 



Note, however, that the advantage of considering the concept of recurrences defined 
in terms of metric distances in phase space instead of correlations is that it allows 
for creating networks based even on individual states without any embedding or 
consideration of groups of states. On the one hand, this independence from a particular 
embedding may be beneficial when dynamical invariants of the studied system are 
of interest. On the other hand, the statistical properties of the resulting recurrence 
networks refiect exclusively the invariant density of states in phase space (in terms 
of certain higher-order statistics), because time-ordering information is lost in this 
framework. Hence, it is hardly possible to distinguish between, e.g., deterministic 
and stochastic dynamics. Here, additional embedding might in fact provide a feasible 
solution to the corresponding identification problem. 

§ As an alternative measure, one could consider the maximum distance of the fc-th nearest neighbour 
as a measure for phase space density. 




(4) 
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Following the above considerations, it can be argued that the concept of recurrence 
networks yields a general framework for transferring time series into complex networks 
in a dynamically meaningful way. In particular, this approach can be applied (i) to 
both univariate as well as multivariate time series (phase space trajectories) (ii) with and 
without pronounced oscillatory components and (iii) with as well as without embedding. 
Consequently, unlike for most existing techniques, there are no fundamental restrictions 
with respect to its practical applicability to arbitrary time series. 

While the definitions of edges and vertices in our approach have already been given 
above (Table [1]), we now provide a geometrical interpretation of a third important 
network entity, the path, within the framework of recurrence networks (Table [2]). A 
path between two vertices i to j in a simple graph without multiple edges can be 
written as an ordered sequence of the vertices it contains, i.e., (i, fci, . . . , /cj- ^._i, j), 
where the associated number of edges measures the length of the path. In phase 
space, a path in the recurrence network is hence defined as a sequence of mutually 
overlapping e-balls i?e(xi), i?e(xfcj, . . . , i?,(xfc;_ ._ J, i?,(xj), where -B,(xi) n i?,(xtfcj 7^ 
0,...,i?,(x, jni?,(x,.)^0i. 

Due to the natural interpretation of vertices, edges and paths, the topological 
characteristics of a recurrence network closely capture the fundamental phase space 
properties of the dynamical system that has generated the considered time series. In 
the following, we will present a detailed analysis of the corresponding analogies for 
different network properties that are defined on a local (i.e. considering only the direct 
neighbourhood of a vertex), intermediate (i.e. considering the neighbourhood of the 
neighbours of a vertex), and global (i.e. considering all vertices) scale (Table [3])0. It 
has to be emphasised that these quantities can be considered as (partly novel and 
complementary) measures in the framework of RQA. 

3.1. Local network properties 

3.1.1. Degree centrality (local recurrence rate). As a first measure that allows to 
quantify the importance of a vertex in a complex network, the degree centrality of 
a vertex v, k^, is defined as the number of neighbours, i.e. the number of vertices i v 
that are directly connected with v. 

N 

ky = ^ ] Ay^j. (5) 

i=l 

Note that in general, the sum is taken over all i ^ v. However, according to our definition 
([2]), we skip the corresponding condition in the following. Normalising this measure by 

II An e-ball centered at state vector x is defined as the open set i?e(x) = {y € R™ : ||x — y|| < e}. 
^ Alternatively, one may classify the corresponding phase space properties according to the fact whether 
they refer to individual points, small regions, or the entire phase space. In this respect, measures related 
to a single vertex v (centralities, local clustering coefficient, local degree anomaly) give local, those 
related to a specific edge (i, j) (shortest path length, matching index, edge betweenness) intermediate, 
and all others global information about the phase space properties. 
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Scale Recurrence network Phase space 



Local Edge density p 

Degree centrality 

Intermediate Clustering coefficient C 

Local degree anomaly Ak, 
Assortativity TZ 
Matching index /i^j 

Global Average path length C 

Network diameter V 
Closeness centrality Cy 
Betweenness centrality by 



Global recurrence rate RR 
Local recurrence rate RRy 
Invariant objects 
Local heterogeneity of phase space density 
Continuity of phase space density 
Twinness of i,j 
Mean phase space separation {dij).j 
Phase space diameter A 
Local centeredness in phase space 
Local attractor fractionation 



Table 3. Correspondence between recurrence network measures and phase space 
properties. Specific terms are discussed in the text. 



the maximum number of possible connections, N — 1, one gets the local connectivity 



1 ^ 

— Y^Ay,, = RRy, (6) 



which, from the recurrence plot point of view, corresponds to the local recurrence rate 
RRy of the state v. Thus, the degree centrahty and local connectivity yield an estimator 
for the local phase space density, since for a vertex v located at position in phase 
space, 

hky{e) + 1)^ f dx p(x) ^ (2e)Xx,) (7) 

(when using the maximum norm) and, hence, 

p(x.)=linijim A±l. (8) 

In complex network studies, one is often interested in the frequency distribution of 
degree centralities, P{k), in particular, the presence of an algebraic scaling behaviour, 
which is characteristic for scale-free networks [20]. However, although several authors 
have recently focussed their attention on this characteristic obtained from different types 
of complex networks derived from time series [30, 34-36, 39, 44-47, 54-56] , we would hke 
to underline that for a complete characterisation of the phase space properties of a 
dynamical system, one should prefer studying not only degree centralities, but also 
other higher-order statistical measures. 

3.1.2. Edge density (global recurrence rate). In some situations, it is useful not to 
consider the full distribution of degree centralities in a network, but to focus on the 
mean degree of all vertices 
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as a simple characteristic quantity of this distribution, where 




N{N -1) 



(10) 



■2 



is the total number of edges in the recurrence network. The mean degree centrality 
{k) is directly proportional to the edge density p of the network or, alternatively, its 
recurrence plot equivalent, the global recurrence rate RR, 



Note that the recurrence rate coincides with the definition of the correlation integral 
C2(e), which is commonly used to estimate the correlation dimension D2, for example, 
using the Grassberger-Procaccia algorithm [59]. 

The connection between the edge density and the correlation dimension can be 
understood by the fact that the local recurrence rate RR^ of a vertex v corresponds 
to the measure of a m-dimensional ball -Be(x^) of radius e centered at the point x^, 
in the m-dimensional phase space in the limit that time goes to infinity {N — > 00). 
When considering the Euclidean norm as a distance measure in phase space, these 
balls are defined as hyperspheres, for the maximum norm as hypercubes etc. Then, 
the pointwise (information) dimension of the probability measure p at is defined 
as Dpl'Xy) = —lim^^Q {In fi{B^{:Kv))/ \ne) [A]. Due to the heterogeneity of the phase 
space visited by the trajectory (i.e., the non-uniform phase space density that results 
in different degree centralities k^ in different parts of this space), the proper estimation 
of Dp is a nontrivial task and often requires expensive computational power and a 
high data quality and quantity. Thus, one may expect a better statistics for D2, since 
it more heavily weights regions of the phase space which have a higher probability 
measure fi. Though the correlation integral has been well established in the literature 
for estimating the correlation dimension, we point out the improvement in estimating 
D2 based on the diagonal lines in Rij{e), which yields an algorithm that is independent 
of the embedding parameters [17]. Consequently, the recurrence network representation 
Aij of a time series fully conserves the geometric properties of the phase space of the 
underlying dynamical system. 

3.2. Intermediate scale network properties 

3.2.1. Local clustering coefficient. The clustering coefficient of a vertex Cy, 
characterises the density of connections in the direct neighbourhood of this vertex in 
terms of the density of connections between all vertices that are incident with v. In 




(11) 



RR{e) = C2(e) ~ e 



-D2 
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many networks, such loop structures formed by three vertices occur more often than 
one would expect for a completely random network. Hence, high clustering coefficients 
reveal a specific type of structure in a network, which is related to the cliquishness of a 
vertex [22]. 

In this work, we consider the definition of clustering coefficient proposed by Watts 
and Strogatz [22], 

C. = TTTi:^^.'', (12) 

where is the total number of closed triangles including vertex v, which is bound by 
the maximum possible value of ky{ky — l)/2. For vertices of degree = or 1 (isolated 
or tree-like points, respectively), the clustering coefficient is defined as Cy = 0, as such 
vertices cannot participate in triangles by definition. 

Equation ( |T2i) can be rewritten in terms of conditional probabilities as 

Cv = P{^i,j = i| = 1, Ay J = 1) = — ^ 2 ~r\ ^^"^^ 

using Bayes' theorem, with 

I.. . = 1^ = 

{N-1){N -2 



^ N N 

P{Ay^i = l,AyJ = 1) = 1 V/V - 9^ 5Z X] (14) 



i=l j=l,j^i 



and a similar expression for P{Aij = 1, Ay^i = 1, Ayj = 1). As for a recurrence network, 
the value of Aij depends only on the phase space distance and the choice of e, the latter 
relationship may be used to derive analytical results at least for one-dimensional systems 
based on their invariant density. Corresponding details can be found in Appendix A 



including the corresponding treatment of the Bernoulli and logistic maps as specific 
examples. 



3.2.2. Global clustering coefficient. As for the degree centrality, one may consider the 
average value of the clustering coefficient taken over all vertices of a network, the so- 
called global clustering coefficient 

1 ^ 

as a global characteristic parameter of the topology of a network. One expects that the 
value of C is - for a given dynamical system with a phase space density p(x) - in the 
asymptotic limit N —>■ 00 exclusively determined by the choice of e, which determines 
the scale of resolution. A more detailed discussion of the corresponding effects and their 
implications for certain model systems will be given in Section HI 
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3.2.3. Mean nearest neighbour degree. The mean nearest neighbour degree A;"" of vertex 
V gives the average degree in the neighbourhood of v, 

1 ^ 

^r^TrE^''.^^- (16) 

i=i 

The degree centrahty k^ is a measure of the density of states in the immediate 
neighbourhood of state v, whereas A;"" can be interpreted to indicate the mean density of 
states in the next neighbourhood (next topological shell of neighbours) of state v. Hence, 
both measures taken together contain information about the local density anomaly in 
the vicinity of v, which we propose to measure by the local degree anomaly 

AA;, = k,- A;^. (17) 

Vertices with a positive degree anomaly (Aky > 0) hence indicate local maxima of phase 
space density, while such which Aky < correspond to local density minima. Hence, 
the local degree anomaly may be considered as a proxy for the local heterogeneity of 
the phase space density. In a similar way, the average absohitc value of the local degree 
anomaly, (|AA;^,|)^, serves as a measure for the overall spatial heterogeneity of the phase 
space density profile. 



3.2.4- Assortativity. A network is called assortative if vertices tend to connect 
preferentially to vertices of a similar degree A;. On the other hand, it is called 
disassortative if vertices of high degree prefer to connect to vertices of low degree, and 
vice versa. Hence, assortativity can be quantified by the Pearson correlation coefficient 
of the vertex degrees on both ends of all edges [22,60], 



7^ 



i ^j>i kikjAij ^j>i |(^« + kj)Aij 



1 2 



2 ■ 



(18) 



If the density of states in phase space hardly varies within an e-ball, the degrees on 
either ends of an edge will tend to be similar and hence the assortativity coefficient TZ 
will be positive. This means, that the more continuous and slowly changing the density 
of states is, the closer TZ will be to its maximum value one. Within the framework of 
recurrence networks, TZ can hence be interpreted as a measure of the continuity of the 
density of states or put differently, of the fragmentation of the attractor. Note that this 
aspect has not yet been specifically addressed by other nonlinear measures, in particular, 
within the RQA framework. 



3.2.5. Matching index (twinness). The overlap of the neighbourhood spaces of two 
vertices i,j is measured by the matching index 

ki + kj — Ylii=i 
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where Hij = if there are no common neighbours, and /ijj- = 1 if the neighbourhoods 
coincide [22]. Using the notion of e-balls around points in phase space, one may 
alternatively write 

. ^ /.(i?e(x,.)ni?.(x,)) 

/i(5e(x.)) + /i(Se(x.)) - /i(5.(x,) n 5,(X,)) ' ^ ' 

Due to the spatial constraints of the recurrence network, the neighbourhood spaces 
of i,j can only overlap if 

dij = ||xj — Xjll < 2e, (21) 

i.e., fiij = for all pairs of vertices with dij > 2e. Moreover, the matching 

index /Zjj decreases on average with an increasing spatial distance dij between the two 
considered states. Note that since Aij = already for dij > e, there may be unconnected 
points with a matching index fiij > 0. 

The matching index of pairs of vertices in a recurrence network is closely related to 
the concept of twins [61], which has recently been successfully applied for constructing 
surrogate data (twin surrogates) in the context of statistical hypothesis testing for the 
presence of complex synchronisation [62,63]. Twins are defined as two states of a 
complex system that share the same neighbourhood in phase space, i.e., the two vertices 
of the recurrence network representing these states have a matching index fiij = 1. 
Hence, the matching index can be used for identifying candidates for twins. Note that 
pairs of vertices i and j in a recurrence network with fii j < 1 can still be considered as 
potential twins, since /ij j = 1 may in some cases be approached by only slight changes of 
the threshold e. Consequently, we suggest interpreting the matching index as a measure 
of the twinness of i and j. Furthermore, it should be noted that adjacent pairs of edges 
{i,j) {Aij = 1) with a low matching index /ijj ~ connect two distinct regions of the 
attractor and may therefore be indicative of geometrical bottlenecks in the dynamics 
(cf. our discussion of the betweenness centrality in Section [3.3.5p . 



3.3. Global network properties 

3.3.1. Shortest path length. As we consider recurrence networks as undirected and 
unweighted, we assume all edges to be of unit length in terms of graph (geodesic) 
distance. Consequently, the distance between any two vertices of the network is defined 
as the length of the shortest path between them. Note that time information is lost after 
transforming the trajectory into a network presentation. Therefore, the terminology of 
the shortest path length /j .,• in the recurrence network reflects the minimum number of 
edges that have to be passed on a graph between a vertex i to a vertex j. In the same 
spirit, /jj is related to the distance of states i and j in phase space. 

In order to better understand the meaning of shortest path lengths, let us study 
their calculation for two toy model series: First, we consider a periodic trajectory 
x{t) = sin(67r ■ O.lt), y{t) = cos(67r ■ O.lt), with t = 0, 1, ■ ■ ■ , 10, i.e. there are = 10 
points in the phase space (Fig. [D^a)). The corresponding recurrence plot is shown for 
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Figure 1. Schematic representation of the transformation of (a) a periodic 
trajectory in phase space (arrows indicate the temporal order of observations) 
into (c) the recurrence network (hnes illustrate the mutual neighbourhood 
relations), (b) gives the associated recurrence plot representation. Table (d) 
lists the resulting shortest path lengths between any two vertices. 



e = I (Fig. [11(b)). As it has already been mentioned above, the recurrence matrix 
Rij and the adjacency matrix Aij of the associated recurrence network are basically 
equivalent. Adopting a common visualisation of connectivity patterns from the literature 
on complex networks, we illustrate the recurrences of the considered model time series by 
placing the individual observations (vertices) on a circle with equal common distances 
(Fig. [11(c)). In this representation, the shortest path length (in the network sense) 
between two vertices i and j corresponds to the smallest number of "jumps" in phase 
space via pairs of neighbours (i.e. recurrences) in phase space. Obviously, the number of 
such jumps is determined by the prescribed value of e and the spatial distance between 
i and j. For instance, the shortest path from vertex z = 1 to j = 10 is h^io = 3 as 
indicated by the matrix of mutual shortest path lengths (Fig. [IJd)). Note that this 
list is symmetric by definition, i.e. lij = Ij^i). The same heuristic analysis can also be 
performed for a general nonperiodic trajectory in phase space as shown in Fig. [2l(a-d). 
We wish to underline that the terminology of shortest path lengths in networks 
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Figure 2. Same as in Fig. [T] for a general non-periodic trajectory 
(3,2,4,8,2,5,7,6,0,2) embedded in a two-dimensional phase space {X,Y) = 
{xt,xt+i). 



does not have a direct relevance to the dynamical evolution of the observed system. 
In contrast, kj measures distances in phase space (among a discrete set of points on 
the attractor) in units of the neighbourhood size e. For example, in the periodic case 
displayed in Fig. [H it takes 9 iterations (time points) from vertex 1 to 10 in the time 
domain, while the shortest path to cover the phase space distance has only a length 
of /i,io = 3. Hence, shortest paths do not allow to infer the temporal evolution of the 
system. Even more, for the path concept in a recurrence network, no information about 
the temporal order of the individual observations is considered (for example, the shortest 
path between vertices 1 and 2 in Fig. [T]is given by the sequence (1,8,5,2), which is not 
ordered in time). 

One should note that if the phase space is strongly fragmented (for instance, in 
the period-3 window of the logistic map, which has been discussed elsewhere [64], the 
phase space consists of three discrete points), the resulting recurrence networks may be 
composed of different disconnected clusters. Furthermore, there might be more than one 
shortest path connecting two nodes. For example, in the aperiodic example in Fig. [21 
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the shortest path from node 1 to node 7, (/ij = 3 as shown in Fig. [21(d)), can be obtained 
by three different choices, that are, (1,2,6,7), (1,5,3,7), and (1,5,6,7). 

3.3.2. Average path length. The average path length C is defined as the mean value of 
the shortest path lengths U^j taken over all pairs of vertices (i, j), 

^=(^-) = Mi73I)E^-- (22) 

Here, for disconnected pairs of vertices, the shortest path length is set to zero by 
definition. Note that in most practical applications, this has no major impact on the 
corresponding statistics. 

The average phase space separation of states {dij) serves as an e-lower bound to £, 
since 

dij < ehj (23) 
due to the triangular inequality, and hence 

{d,j) < eC. (24) 

Interpreted geometrically, this inequality holds because C approximates the average 
distance of states along geodesies on the recurrence network graph (which can be 
considered as the geometric backbone of the attractor) in multiples of e, while (dij) 
gives the mean distance of states in M"* as measured by the norm || ■ ||. 



3.3.3. Network diameter. By a similar argument as used in Eq. (1231) for the average 
path length, the diameter 

T> = max/jj (25) 

of the recurrence network (i.e. the maximum path length) serves as an e-upper bound 
to the estimated diameter 

A = max di J (26) 

«j 

of the attractor in phase space: 

A < eV. (27) 

3.3.4. Closeness centrality. The inverse average shortest path length of vertex v to all 
others in the recurrence network is measured by the closeness centrality [65] 

A^- 1 

c. = ;v , • (28) 

If i and j are not connected, i.e., Aij = 0, the maximum shortest path length in the 
graph, — 1, is used in the sum by definition. In a recurrence network, can be 
geometrically interpreted as measuring the closeness of v to all other states with respect 
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to the average length (in units of e) of geodesic connections on the recurrence network 
graph. In other words, is large if most of the other vertices are reachable in a small 
number of e-jumps from state to state. 

From Eqs. (!23|) and (l28l) . we can see that the inverse closeness is bounded 
from below by the average phase space distance of vertex (state) v to all other vertices 
(states) {dy^i)^ in units of e (geometrical closeness), as measured by the norm || ■ ||, 

;^K.).<c.-^- (29) 
Put differently, geometrical closeness provides an upper bound for topological closeness, 

Cv < e ^ {dv,i)i ■ (30) 

3.3.5. Betweenness centrality. The betweenness centrality 6„ has been originally 
introduced for characterising the importance of individual vertices for the transport 
of information or matter in general complex networks [65] . Unlike the degree centrality 
ky, it is defined locally but depends on global adjacency information. 

Let us assume that information travels through the network on shortest paths. 
There are cjjj shortest paths connecting two nodes i and j. We then regard a node v to 
be an important mediator for the information transport in the network, if it is traversed 
by a large number of all existing shortest paths. Betweenness centrality is given by 

N . , 



E (31) 



tj^v '-' 

where cTjj(f) gives the number of shortest paths from i to j, that include v. Here 
the contribution of shortest paths is weighted by their respective multiplicity aij, the 
physical rational for this normalisation being that the total volume of information flow 
between two vertices, when summed over all shortest paths connecting them, should 
be the same for all pairs in the network. Hence, in addition to degree and closeness 
centralities, betweenness centrality yields another possibihty to identify especially 
relevant vertices. 

For a recurrence network, the notion of information transfer is not useful anymore. 
However, one may still argue in a geometric way that high betweenness states are 
typical for regions of sparse phase space density that separate different high-density 
clusters (refering to the information flow analogy mentioned above, one may consider 
the corresponding vertices as geometric bottlenecks). Thus, the occurrence of high 
betweenness values can be a sign of highly fractionated attractors (on the scale resolved 
by the considered threshold e). A more detailed discussion of the corresponding 
implications for some simple model systems will be given in Section HI 



3.3.6. Edge betweenness. While betweenness centrality refers to vertex properties of a 
network, one may define an equivalent measure also based on the number of shortest 
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paths on the network that include a specific edge We refer to the corresponding 

property as the edge betweenness bij. Note that though there is a conceptual difference 
between vertex-related and edge-related betweenness, both quantities are indicators for 
regions of low phase space density that separate regions with higher density (or, to say 
it differently, of regions of high attractor fractionation) and thus have practically the 
same dynamical meaning. 

4. Examples 

In the following, we will show the potentials of the network-theoretic measures discussed 
in the previous section for recurrence networks obtained from three paradigmatic chaotic 
model systems. 

4.I. Model systems 

Basic results for one-dimensional maps have already been described for the logistic 
map (see [64]) based on numerical calculations and are supplemented by some further 
computations in the appendix. At this point, we prefer to discuss in some more detail 
the properties of systems that are defined in somewhat higher dimensions. In particular, 
we consider the Henon map 

Xi+i =yi + l- 1.4x- , yi+i = 0.3xi (32) 

as an example for a chaotic two-dimensional map, and the Rossler system 

^ (x, y, z) = [-y-z,x + 0.2y, 0.2 + z{x - 5.7)) (33) 
as well as the Lorenz system 

^ (x, y, z) = (^10{y - x),x{28 - z) -y,xy- ^z^ (34) 

as two examples for three-dimensional chaotic oscillators. In all following considerations, 
no additional embedding will be used. Note, however, that for the continuous systems, 
temporal correlations between subsequent observations have been excluded by removing 
all sojourn points [66]. 

Figs. [3] and m show examples of typical trajectories of these three model systems. In 
addition, the shortest paths between the first and last point of the individual realisation 
are indicated, underlining the deep conceptual differences between the concepts of 
trajectory (in phase space) and path (in a recurrence network, see Section [3]). 

4-2. e-dependence of global network measures 

Let us first consider the dependence of the global network measures £, C and 71 on the 
choice of the threshold e for our three model systems. 



Recurrence networks - A novel paradigm for nonlinear time series analysis 



21 




Figure 3. One trajectory of the Henon map (attractor indicated by grey 
dots) formed by 100 iterations (i.e. N = 101), indicated by circles. The 
initial condition is marked by a square the size of which corresponds to 
the neighbourhood threshold e = 0.25 (maximum norm) considered in the 
derivation of the corresponding recurrence network. The shortest path (with 
^1,101 = 10) between the first and the last point is indicated by a continuous 
red line. 



The variations of C with the threshold e are shown in Fig. O and verify the existence 
of an inverse relationship of a corresponding lower bound postulated in Eq. (^^. 

For the global clustering coefficient, the dependence on e is more complicated and 
depends on the specific properties of the considered system (Fig. [5]). In particular, 
while for too small e, problems may occur, since the recurrence network may decompose 
into different disconnected clusters for a length of the considered time series, for 
intermediate threshold values, an approximately linear increase of C with e seems to be 
a common feature of all three examples. Following the discussion of the behaviour of 
one-dimensional maps in Appendix A we may argue that this increase is most likely 
related to the effect of the attractor boundaries. 

Finally, concerning the assortativity coefficient TZ, we observe that for small e, 
the recurrence networks are highly assortative (e.g. TZ is close to 1). This behaviour 
can be related to the fact that in case of small neighbourhoods, these phase space 
regions are usually characterised by only weak variations of the phase space density, 
so that neighbouring vertices have a tendency to obey a similar degree. As e becomes 
larger, larger regions of the phase space are covered, where the density may vary much 
stronger, which implies that the degrees of neighbouring vertices become less similar. 
Note, however, that since in this case, the mutual overlap of the different neighbourhoods 
becomes successively larger, there is still a significantly positive correlation between the 
degrees of neighbouring vertices. One may further observe that the decrease of TZ with 
e may be interrupted by intermediate increases, which are probably related to some 
preferred spatial scale of the separation of certain dynamically invariant objects such as 
unstable periodic orbits (UPOs). We will come back to this point in Sec. 14.3.31 



Recurrence networks - A novel paradigm for nonlinear time series analysis 22 




20 -20 -20 -10 10 20 

X 



Figure 4. (a,b) Chaotic attractor of the Rossler system (dark dots) and 
reahsation of one particular non-periodic trajectory (blue line), corresponding 
to T = NAt = 11.4 time steps. The size of the considered neighbourhood e = 1 
(maximum norm) is indicated by a red square around the initial condition. 
(Note that the square gives an idealised representation of the considered 
neighbourhood.) The red line displays the shortest path between the initial 
condition and the final value on this trajectory (I = 5). (c,d) One example 
trajectory of the Lorenz system with T = 5 time steps, and resulting shortest 
path between initial and final state (e = 1.5, Z = 9). 

4-3. Spatial distributions of vertex properties 

In the following, we will study the interrelationships between local network properties 
and structural features of the phase space for the three considered chaotic model systems. 

4-3.1. Degree centrality. When considering the degree centrality or, equivalently, 
the local density for all vertices of the network, a broad range of variability is found 
(Fig. E]). In particular, the behaviour follows the expectation that regions with a high 
phase space density (for example, the merger of the two scrolls of the Lorenz oscillator) 
also reveal a high density of vertices and, hence, high degree centralities. Note that the 
calculation of a recurrence plot depends on the parameter e, which should be tailored 
to the considered system under study and the specific questions one wishes to address. 
Several "rules of thumb" for the choice of the threshold e have been advocated in the 
literature [9, 10]. It has been suggested that the choice of e to achieve a fixed recurrence 
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Figure 5. Dependence of (a), (b), (c) the average path length £, (d), (e), (f) 
the global clustering coefficient C, and (g), (h), (i) the assortativity coefficient 
TZ for the Henon map, Rossler, and Lorenz system (from left to right) using the 
maximum norm. Dashed lines in the plots on £(e) indicate the approximate 
presence of the theoretically expected 1/e dependence of the average path 
length. Note that although a normalised threshold e might yield a better 
comparability of the results for the different systems, we prefer using the 
absolute values here since the typical normalisation factors - either an empirical 
estimate of the standard deviation of the phase space density (which may be 
crucially influenced by asymmetric densities) or the attractor diameter (which 
is itself not a priori known in advance) have certain conceptual problems. 

rate RR is helpful for the estimation of dynamical invariants in many systems [9]. 
Therefore, this procedure will be adopted here to obtain an overall visualisation of 
the degree centrality in phase space, with RR = p ^ 0.03 (which lies within the 
typical scaling region of the correlation integral). However, as we will see later, for the 
local clustering coefficients (Sec. I4.3.3P disclosing local fine structures of the phase space 
density, it is necessary to choose smaller e. 

4-3.2. Closeness centrality. Figure [7] refiects the spatial distribution of the closeness 
centrality Cy . In good agreement with our previous theoretical considerations on the 
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Figure 6. Colour-coded representation of the local recurrence rate RR^ 
(proportional to the degree centrality ky) in phase space (a) Henon map (A^ = 
10,000), (b) Rossler {N = 10,000), and (c) Lorenz systems {N = 20,000). The 
value of e for each case is chosen in such a way that the global recurrence rate 
RR = 0.03. 
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Figure 7. Colour-coded representation of the closeness centrality (Eq. (p8]) ) 
in phase space (a) Henon map, (b) Rossler, and (c) Lorenz systems {N and e 
as in Fig. [6|). 



geometric meaning of this measure (Sec. 13.3.41) . we find high values of near the centre 
of gravity of the attractor in phase space, and low values at phase space regions that 
have large distances from this centre. 

4-3.3. Clustering coefficient. Concerning the local clustering coefficient, one may 
suppose that in the case of high-density regions in phase space, there are many vertices 
located in the vicinity of a specific vertex, in particular, in a distance that does not 
exceed e/2. By definition, these vertices must then also be adjacent to each other, 
which gives considerable contributions to the clustering coefficient. In contrast, for low- 
density regions, one may argue that even if there are more than one vertices in some 
e-neighbourhood of a vertex (i.e., >2), it is less likely that these are also separated by 
a distance that is smaller than e. Following these considerations, one might expect some 
relationship between the degree centrality and the local clustering coefficient. However, 
as Figs. [HI El and [in] demonstrate, does clearly reveal more and different structural 
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Figure 8. Scatter diagrams between the local recurrence rate RRi, and the 
local clustering coefficient Cv for typical realisations of (a) Henon map, (b) 
Rossler, and (c) Lorenz system. The inserted values ps give the corresponding 
rank-order correlation coefficient (Spearman's Rho). 




Figure 9. (a) Colour-coded representation of the local clustering coefficients 
for the Henon map in its phase space, (b) Relationship between the local 
clustering coefficients and the stable manifold of the Henon map. The finite 
length segment of the stable manifold is calculated by the method described 
in [67,68] (with 20,000 iterations). 

properties than the degree centrality alone (note that here, a smaller value of e has been 
chosen to disclose the local fine structures of the phase space density) . In particular, a 
visual comparison with Fig. reveals that the clustering coefficient characterises some 
specific higher-order characteristics of the phase space density. 

Beside the effect of the local phase space density, we argue that the local clustering 
coefficient does also depend on the spatial filling (i.e., the homogeneity of the phase 
space density) in the neighbourhood of the considered point. In particular, in the 
case of a two-dimensional system, an alignment of vertices along a one- dimensional 
subspace will produce a clearly lower clustering coefficient than a homogeneous filling 
of the neighbourhood. This behaviour is underlined in Fig. [9] for the Henon map, where 
maximum values of = 1 can be particularly found at the two tips of the attractor. 
Hence, the local clustering coefficient can be considered as an entropy-like characteristic 
in that it quantifies the homogeneity of the phase space density in the neighbourhood 



Recurrence networks - A novel paradigm for nonlinear time series analysis 26 




10 



20 



20 -20 20 -20 



Figure 10. (a) Colour-coded representation of the local clustering coefficients 
Cy for the Rossler system (A^ = 10, 000) in its phase space, (b) Several locations 
of the unstable periodic orbits with low periods, obtained using a method 
based on the windows of parallel lines in the corresponding recurrence plot 
as described in [9]. The chosen value of e corresponds to a global recurrence 
rate of RR = 0.01. (c),(d) Same as (a),(b) for the Lorenz system (A'' = 20, 000). 

of a vertex. From the theory of spatial random graphs [69], which may be assumed to 
yield the lowest possible clustering coefficients among all spatial networks, it is known 
that for a given dimension m of the considered system, in the asymptotic limit N ^ oo 
and e — i> 0, the possible values of C„ are bound between the corresponding theoretical 
value and 1. Note that this lower bound systematically decreases with increasing m, 
which appears to be reasonable if one interpretes as an entropy-like quantity. To be 
more specific, according to Dall and Christensen [69], this decay follows an exponential 
function for sufficiently large embedding dimensions. 

The presence of distinct structures in the spatial profile of the local clustering 
coefficient is related to the emergence of specific dynamically invariant objects in the 
considered model systems. In the case of Henon map (Fig. [9]), there is a clear tendency 
that points that are close to the stable manifold associated with the system have 
remarkably higher values of Cy. Note, however, that because of finite size effects. 
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this coincidence cannot be found for all corresponding regions of the phase space. For 
the two continuous systems (Fig. [TOl) . points close to the trapping regimes of UPOs 
have higher clustering coefficients. It is, in some sense, trivial to understand the role 
of UPOs in forming such regimes of higher clustering. Whenever the trajectory of 
the corresponding systems visits the neighbourhood of an UPO, it is captured in this 
neighbourhood for a certain finite time, during which the probability of recurrences 
is increased. Furthermore, once the trajectory is trapped, the local divergence rate 
becomes smaller. This smaller local divergence rate is captured by the clustering 
coefficient (in terms of higher-order correlations between neighbours of a vertex) . As for 
the finite-e effect in the Henon map, the regions with increased clustering coefficients in 
most cases only coincide with UPOs of lower periods. Therefore, in Fig. [TOl only a few 
UPOs of low order are shown for comparison. Note that if two UPOs are separated by a 
distance smaller than e in phase space, the clustering coefficient is not able to distinguish 
between these two structures and, hence, shows a broad band with increased values. 
Following this argumentation, C„ is a useful measure for detecting phase space regions 
with a high density of low-order UPOs, which is in good agreement with corresponding 
considerations in [39]. 

4.3.4- Betweenness centrality. Our interpretation of the betweenness centrality in 
Sec. 13.3.51 implies that is a rather sensitive measure of the local fragmentation of 
the attractor and thus may give complementary information especially on very small 
scales. Unfortunately, numerical limitations in the calculation of this measure did not 
allow us to explore the limit of small neighbourhoods (e —>■ 0). However, from our 
computations with somewhat larger thresholds (see Fig. [TTj) . we can already derive some 
general statements about the behaviour of betweenness centrality for the considered 
model systems. First, note that regions close to the outer boundaries of the attractor 
(in contrast to those in the vicinity of the inner boundaries, e.g., of the Rossler oscillator) 
are not important for many shortest path connections on the recurrence network. Hence, 
vertices settled in the corresponding parts of the phase space are characterised by low 
betweenness values. Second, if there are pronounced regions with rather few isolated 
points in between high-density regions (for example, between two UPOs in the Rossler or 
Lorenz systems), there is an increasing number of shortest paths crossing these vertices, 
which leads to higher values of b^. In turn, vertices in the vicinity of UPOs (i.e., 
high-density regions) show lower betweenness values. Therefore, betweenness centrality 
provides a complementary view on the attractor geometry in comparison to the local 
clustering coefficient (Fig. [TOl). 

4.4- Spatial distributions of edge properties 

Similarly to the local vertex properties, one may also study the characteristics of different 
edges in the recurrence networks. Since the resulting structures are more pronounced 
than for the three model systems considered so far. Fig. [12] shows the matching index 
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Figure 11. Logarithm of the betweenness centrahty for (a) Rossler and (b) 
Lorenz system (A^ and e as in Fig. llU|) . Points shown as circles have betweenness 
values below the lower limit of the displayed colour scale. 

and edge betweenness for one realisation of the logistic map Xi+i = axj(l — Xi) in the 
intermittent chaotic regime (see [64]). The presence of intermittent dynamics can be 
clearly seen from the recurrence plots in terms of extended square recurrence patterns, 
which hence lead to mutually connected vertices of the associated recurrence networks 
that correspond to subsequent points in time. 

Figures [T2] and [T3] show the complex dependence between phase space distance dij, 
matching index /Xj and edge betweenness bij. For the matching index, the results are 
consistent with our theoretical considerations presented in Sec. 13.2.51 In particular, 
for dij 0, we have fiij 1, while for dij 2e, fiij 0. Concerning the 
temporal evolution during the laminar (intermittent) phase, one may recognise that 
at the beginning, there is hardly any change in the state of the system, hence, dij is 
very small for subsequent points in time (vertices of the recurrence network), which 
relates to large values of the matching index near 1. As the laminar phases are close 
to their termination, chaotic variations emerge and rise in amplitude, which leads to a 
subsequent increase of dij and, hence, decrease of /ijj. 

Concerning the edge betweenness bij (the spatial pattern of which is very similar to 
that of the vertex-based betweenness centrality b^ due to the spatial proximity of edge 
and corresponding vertices), the overall behaviour is opposite to that of /ij j. During 
laminar phases, we find that since all states are very close to each other, possible shortest 
connections may alternatively pass through a variety of different edges, leading to low 
values of the edge betweenness. Close to the termination, there is in turn an increase of 
this measure. However, the most interesting feature of the edge betweenness is presented 
by isolated edges with very high values of bij, which correspond to rarely visited phase 
space regions between intervals of higher phase space density. More specifically, the 
average edge betweenness of vertices in such low-density regions may exceed that of 
high-density regions by orders of magnitude (Fig. [T2ll . 
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Figure 12. (a,d) One example trajectory of the logistic map at a = 3.679 
(chaotic regime with intermittency) , e = O-OlSa^;, N = 1000 (only a part of the 
trajectory with 250 points is shown), (b) Matching index ^Ujj- between all pairs 
of vertices (colour-coded), (c) Average matching index /ij of all vertices of the 
considered recurrence network. (e,f) As in (b,c) for the logarithm of the edge 
betweenness bij. 
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Figure 13. Scatter diagram of the matching index fiij against the phase space 
distance dij for the logistic map at a = 3.679 (parameters as in Fig. [T2]) . 
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5. Conclusions 

This paper has reconsidered the analysis of time series from complex systems by means 
of complex network theory. We have argued that most existing approaches for such an 
analysis suffer from certain methodological limitations or a lack of generality in their 
applicability. As an appealing solution, we have suggested recurrence networks as a 
unifying framework for studying time series as complex networks, which is based on 
the idea of recurrence plots. As we have argued, this specific approach is applicable to 
univariate as well as multivariate time series without or with embedding. In addition, 
recurrence networks can be applied for studying time series with non-equidistant time- 
scales and allow the construction of simple significance tests with respect to the 
associated network-theoretic measures [64]. 

As a main achievement, we have provided a thorough reinterpretation of a variety 
of statistical measures from network theory computed for recurrence networks in terms 
of phase space properties of dynamical systems. Since all time ordering information is 
lost in this approach, all complex network characteristics are dynamically invariant, i.e., 
they are only sensitive to certain properties of the invariant density of the considered 
dynamical system. From this invariance, it follows that specific measures such as the 
local clustering coefficient may be used for detecting dynamically invariant objects like 
unstable periodic orbits or chaotic saddles. On the other hand, this feature also implies 
that the proposed method cannot be used to distinguish between deterministic (chaotic) 
and stochastic systems, which is exemplified by our comparison between Bernoulli map 



and uniform noise in [Appendix A[ As a possibility to overcome this potential point 
of criticism to our method, we emphasise that an additional embedding should change 
the properties of deterministic systems in a different way than for a stochastic system, 
hence, studying complex network properties in dependence on the embedding dimension 
might help solving this interpretation problem. We will further elaborate this idea in 
future research. 

Using wide-spread statistical characteristics of complex networks such as the 
"trinity" of centrality measures (degree, closeness, betweenness) and the clustering 
coefficient, we were able to provide a detailed interpretation of the corresponding results 
for recurrence networks in terms of higher-order phase space properties. In particular, 
degree centrality relates to the local density, closeness centrality to the average 
geometrical proximity of an observation to all other observations, and betweenness 
centrality to the local fragmentation of points in phase space. For the clustering 
coefficient, the consideration of different model systems has revealed that apart from 
possible density effects, the clustering properties are related to the spatial filling of the 
phase space, which become important close to the attractor boundaries and certain 
dynamically invariant objects such as invariant manifolds or unstable periodic orbits. 
Additionally, we have presented a rigorous analytical treatment of the local and global 
clustering coefficients of the recurrence networks of one-dimensional chaotic maps. 



which perfectly matches our numerical results (see Appendix A). A possible further 
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relationship between the local clustering coefficient of a recurrence network and the 
local Lyapunov exponent of the underlying dynamical system remains a topic for future 
work. 

With respect to existing recurrence plot based methods of time series analysis, 
e.g., recurrence quantification analysis (RQA), we would like to emphasise that our 
approach yields a complementary view on the phase space properties of the underlying 
dynamical system. In particular, one should note that nearly all of the considered 
network-theoretical measures have no direct equivalents in traditional RQA and vice 
versa. Hence, one may think of situations when either one of the two frameworks 
(i.e., RQA or recurrence networks) may provide better results than the other. In turn, 
there have been some recent approaches on applying methods of time series analysis 
to general complex networks (e.g. [70,71]). We would like to underline that due to 
the duality between the recurrence matrix of a time series and the adjacency matrix of 
the associated recurrence network, RQA might be another promising candidate for this 
purpose (as long as one restricts oneself to measures that are invariant under re-ordering 
of time), which might yield interesting complementary insights on complex networks in a 
variety of different situations. A more detailed investigation of a corresponding approach 
will be subject of future research. 
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Ministry for Science and Education. 



Appendix A. Clustering coefficient of recurrence networks for 
one-dimensional maps 

In order to compute the clustering coefficient of the recurrence network of a dynamical 
system, certain system-specific integrals have to be solved. For the case of one- 
dimensional maps that are defined on the integral [0, 1], these integrals can be explicitly 
expressed and eventually also evaluated analytically. 

Appendix A. 1. General treatment 

The computation of the local clustering coefficient involves the computation of certain 
integrals to estimate the required probability terms discussed in Section |3l In particular, 
we have that for a given vertex f at a point x^, 

-Pi^Ay^i 1: -^v,j 1) -^(1 1'^v '^i II ^ ^; I \'^v '^j 1 1 ^ ^) 

= P{\\xy -Xi\ \ < e)P(||a;^ - Xj\ \ < e) 



2 



dx p{x) 



(A.l) 



P(\\x„ - XiW < e) 

.J Xv-e 

for the condition in the denumerator of Eq. (fT3l) . where we have made use of the fact 
that the observations Xi and xj can be considered to be independent of each other. 
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However, this independence condition is not fulfilled in the corresponding numerator 
P{Aij = 1, Am = 1, Ayj = 1), which therefore requires a subtle choice of the integration 
boundaries to assure a correct treatment of the three-point relationships 

rxv+t i'ram{x+t,Xv+e) 

P(Aij = 1, A^^i = 1, A^j = 1) = / dx p{x) / dy p{y). (A.2) 

Jxy—t J ma3i{x—e,Xv—t) 

As one may see in Fig. lAll the correct integration is not trivial. In particular, if we use 
the abbreviations 

Ii{a,b) = / dx p{x) (A. 3) 

J a 

(A.4) 



fb r fd 

l2{a,b;c,d) = / dx p{x) / dy p{y) 

J a L J c 



we obtain the following expressions for the local clustering coefficient Cy = C{xy) if 
e < 0.5: 

< Xy < e: 

_ 12(0, Xy] 0, X + e) + hixy, e; 0, Xy + e) + hie, Xy + e; x - e, Xy + e) 
- h{0,xy + ey 

.(1) .(2) .(3) l^-^J 

_ -'2 "T -'2 "T -'2 

~ /i(0,x, + e)2 ' 
e < Xy < 1 — e: 

h{xv - e, Xy] Xy - e,x + e) + hixy, Xy + e;x - e,Xy + e) 
Ii{xy -e,Xy + e)2 



- e,Xy + e)^' 

1 — e < a;^, < 1: 

_ /2(xt; - e,l - e]Xy - e,x + e) + 12(1 - e, x^,; - e, 1) + /2(x^„ 1; x - e, 1) 



Ji(xt, - e, 1 



\2 



/i(x„-e, 1)2 

Understanding the global clustering coefficient as the expectation value of the local 
one (taken over the whole possible range in x) , one may use the following expression for 
deriving its value: 

C(e) = I dxy p{xy) Cy{xy; e) = I dxy p{xy) ^^^'L . . ^'^''iN — 
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Figure Al. Schematic representation of the integration boundaries entering 
the three-point relationships in Eq. ()A.2p for e < 0.5. 



Appendix A. 2. Bernoulli map 

Among all nonlinear maps defined on the unit interval [0,1], the Bernoulli map 
Xi+i = 2xi mod 1 has the simplest possible invariant density p{x) = 10- This allows 
an easy evaluation of the integrals for an analytic computation of the local clustering 
coefficient, being aware of the restricted integration range. As a result, one finds: 



Cy{x,;e) = < f, e < X < 1 - e . (A.9) 




In particular, for e ^ (and 00), we have C^, — > | V x g]0, 1[, which corresponds to 
the value of random geometric graphs in one dimension [69] . Hence, one may speculate 
about this value being a universal limit for the recurrence networks of one- dimensional 
chaotic maps. Moreover, for x — > and x 1, we have —>■ 1 independent of e. This 
behaviour is consistent with our previous observations concerning the effect of sharp 
attractor boundaries on the local clustering coefficient, for example, in case of the tips 
of the Henon attractor (Sec. 14.3.31) . 

For the global clustering coefficient, the computation of Eq. (lA.SP leads to 



C(e) = ^ + e f41n2- I ) , (A.IO) 



4 V 2 

Note that from the recurrence network properties, this map cannot be distinguished from a stochastic 
process with a uniform density in the same interval (see Fig. IA2p . 
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Figure A2. Dependence of the local clustering coefficient on the threshold 
e for (a) the Bernoulli map and (b) uniformly distributed noise with values 
in [0,1]. Note that due to the same distribution of the data, i.e., the same 
density in the one-dimensional phase space, both curves are equal and match 
the theoretical expectations (dashed red lines). 



which is in good agreement with numerical results (see Figure R2|) . Hence, deviations 
from the theoretical value | occur exclusively due to boundary effects, and may even 
lead to C = 1 for very large thresholds e. 

Appendix A. 3. Logistic map for a = 4 

The logistic map at a = 4 is known to have the invariant density 

1 1 



p{x) 

TT y/x{l - X) 

With this, we find the following expressions: 



— ^yA^j — 1) — 

I — ;^ arcsin(l — 2xv — 2e), < a;^, < e 

^ {arcsin(l — 2xv + 2e) — arcsin(l — 2xv — 2e)} , e < x^ < 1 — e 

I + ^ arcsin(l — 2xy + 2e), 1 — e < a;^, < 1 



(A.ll) 



(A.12) 
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and 



r(l) 



r(2) 



r(3) 



1 1 

arcsinf 1 

4 27r ^ 



1 



I dx 



arcsin(l — 2x — 2e) 
a/^(1 - 2;) 



1 

2^ 



1 



arcsin(l — 2a;„ — 2e) (arcsin(l — Ix^) — arcsin(l — 2e)) 



— arcsinf 1 

77^ 



2xv — 2e) (arcsin(l — 2xv — 2e) 
2x + 2e) 



arcsin(l — 2e)) 



, arcsinCl 
dx - 



X 



r(4) 



77^ 



arcsin(l — 2xv + 2e) (arcsin(l — 2xv + 2e) — arcsin(l — 2xv)) 



- — / rfx 



arcsin(l — 2x — 2e) 



(A.13) 
(A.14) 

(A.15) 

(A.16) 



x] 



r(5) 



77' 



arcsin(l — 2x^ — 2e) (arcsin(l — 2x^ — 2e) — arcsin(l — 2xy)) 
1 



+ 



77^ 



, arcsinf 1 - 2x + 2e) 
dx - 



(A.17) 



x] 



r(6) 



77^ 



arcsin(l — 2xv + 2e) (arcsin(l — 2x^ + 2e) + arcsin(l — 2e)) 



r(7) 



r(8) 



1 
1 

2^ 
1 



l-e 
1 



arcsin(l — 2x — 2e) 
arcsin(l — 2x.u + 2e) 



, H arcsinf 1 

4 277 



2x,) 



77^ 



'arcsin(l — 2x„) + arcsin(l 
arcsin(l — 2a; + 2e) 



2e)) 



(A.18) 
(A.19) 
(A.20) 



Note that the integrals ,L 



(2) ^(3)^ 
) -'2 > 



and (Jf ^ 



r{7) r(6)N 

5 -'2 5 -'2 ^ 

(4) , r(5) 



can be transformed into each 
other by the transformation x ^ 1 — x, and that I"^' + /^^^ is invariant under the same 
transformation, reflecting the corresponding symmetry of the phase space density p{x). 

Since the remaining integrals can only be solved numerically, we are not able to 
give an explicit equation for the functional dependence of the clustering coefficient on 
both X and e. However, comparing the numerical solution of our analytical results 
with the clustering coefficients of recurrence networks of one realisation of the logistic 
map at a = 4, we find (apart from remaining fluctuations due to the flnite length of 
the considered time series) an excellent agreement (see Fig. IA3p . In particular, for 
the attractor boundaries at x = and x = 1, we have = 1 independent of e 
as for the Bernoulli map. Moreover, we observe almost uniform values of the local 
clustering coefficient close to | within the interval [e, 1 — e], and a systematic tendency 
towards larger values close to the boundaries. Since the measure of the latter intervals 
systematically increases with increasing e, these boundary effects are again responsible 
for the systematic increase in the global clustering coefficient C as e becomes larger. 
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Figure A3, (a) Analytical (dashed red line) and simulation results (circles) 
on the dependence of the local clustering coefficient on the spatial position 
X for the logistic map at o = 4 (A^ = 5,000, e = 0.1). (b) As in (a) for the 
dependence of the global clustering coefficient C on e. 
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